#include "stdafx.h"

#include "../MFCD3D11Test.h"
#include "RayTracingHM.h"

//Using CImg library
#include "../Utility/CImg.h"
using namespace cimg_library;

#define  TENDVAL 1000.0f
#define  INITZVALSCALE 0.5f
const float CRayTracingHM::ExtraEpsilon = 2.0e-3f;

CRayTracingHM::CRayTracingHM():m_pScene(NULL), m_iCurHierBldStyle(GreedySplit), m_iNumRaysPerRowIn2D(512*2), m_iVGonioOutputImgWidth(512), 
	m_iVGonioOutputImgHeight(512), m_fTotalIncidentEnergy(1000000.0)
{

}

CRayTracingHM::~CRayTracingHM()
{
	m_arrRays.RemoveAll();
	m_arrTriNormals.RemoveAll();
	SAFE_DELETE(m_pScene);
}

bool CRayTracingHM::createGeometryFrom2DVertArr(const std::vector<VertexNorm>& arrVert, unsigned int iWidth, unsigned int iHeight)
{
	CHiPerfTimer timer;
	timer.Start();

	if(arrVert.size()==0 || (arrVert.size() != iWidth*iHeight))
		return false;

	if(!m_pScene)
		m_pScene = new Scene();
	else
		m_pScene->deleteAll();

	m_arrTriNormals.RemoveAll();

	PoolAllocator *geomPool = m_pScene->getGeometryObjectPool();
	//valid triangle id starts from 1, and 0 is reserved for invalid triangle
	int iID = 1;
	int numTris = 0;
	//prepare triangles for geometry
	for(unsigned int j = 0; j <(iHeight-1); j++)
		for(unsigned int i = 0; i <(iWidth-1); i++)
		{
			unsigned int v1 = j*iWidth+i;
			unsigned int v2 = v1 + 1;
			unsigned int v3 = v1 + iWidth;
			unsigned int v4 = v2 + iWidth;
			float v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z;
			/*********************************************************************/
			//first triangle (v1,v2,v4);
			/*********************************************************************/			
			v0x = arrVert[v1].Position.x; v0y = arrVert[v1].Position.y; v0z = arrVert[v1].Position.z;
			v1x = arrVert[v2].Position.x; v1y = arrVert[v2].Position.y; v1z = arrVert[v2].Position.z;
			v2x = arrVert[v4].Position.x; v2y = arrVert[v4].Position.y; v2z = arrVert[v4].Position.z;

#ifdef DOUBLE_PRECISION_TRIANGLES
			TriangleDP *tri = TriangleDP::poolAlloc(geomPool);
			if(v0z < 0 || v1z < 0 || v2z < 0)
				//invalid triangle with 'No-data' pts
				tri->set(v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z,0);
			else {
				tri->set(v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z,iID++);
				//compute the normal for the triangle
				D3DXVECTOR3 v1(v0x,v0y,v0z);
				D3DXVECTOR3 v2(v1x,v1y,v1z);
				D3DXVECTOR3 v3(v2x,v2y,v2z);
				D3DXVECTOR3 vNorm;
				CD3DUtilityFunction::ComputeTriangleNormal(vNorm,v1,v2,v3);
				m_arrTriNormals.Add(vNorm);
			}

#else
			Triangle *tri = Triangle::poolAlloc(geomPool);
			if(v0z < 0 || v1z < 0 || v2z < 0)
				//invalid triangle with 'No-data' pts	
				tri->set(v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z,0);
			else{
				tri->set(v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z,iID++);
				//compute the normal for the triangle
				D3DXVECTOR3 v1(v0x,v0y,v0z);
				D3DXVECTOR3 v2(v1x,v1y,v1z);
				D3DXVECTOR3 v3(v2x,v2y,v2z);
				D3DXVECTOR3 vNorm;
				CD3DUtilityFunction::ComputeTriangleNormal(vNorm,v1,v2,v3);
				m_arrTriNormals.Add(vNorm);
			}

#endif			
			//add triangle into the scene object			
			m_pScene->addGeometryObject(tri);
			numTris++;
			/*********************************************************************/
			//second triangle (v1, v4, v3)
			/*********************************************************************/
			v1x = arrVert[v3].Position.x; v1y = arrVert[v3].Position.y; v1z = arrVert[v3].Position.z;

#ifdef DOUBLE_PRECISION_TRIANGLES
			TriangleDP *tri1 = TriangleDP::poolAlloc(geomPool);
			if(v0z < 0 || v1z < 0 || v2z < 0)
				//invalid triangle with 'No-data' pts
				tri1->set(v0x,v0y,v0z,v2x,v2y,v2z,v1x,v1y,v1z,0);
			else{
				tri1->set(v0x,v0y,v0z,v2x,v2y,v2z,v1x,v1y,v1z,iID++);
				//compute the normal for the triangle
				D3DXVECTOR3 v1(v0x,v0y,v0z);
				D3DXVECTOR3 v2(v2x,v2y,v2z);
				D3DXVECTOR3 v3(v1x,v1y,v1z);
				D3DXVECTOR3 vNorm;
				CD3DUtilityFunction::ComputeTriangleNormal(vNorm,v1,v2,v3);
				m_arrTriNormals.Add(vNorm);
			}

#else
			Triangle *tri1 = Triangle::poolAlloc(geomPool);
			if(v0z < 0 || v1z < 0 || v2z < 0)
				//invalid triangle with 'No-data' pts	
				tri1->set(v0x,v0y,v0z,v2x,v2y,v2z,v1x,v1y,v1z,0);
			else {
				tri1->set(v0x,v0y,v0z,v2x,v2y,v2z,v1x,v1y,v1z,iID++);				
				//compute the normal for the triangle
				D3DXVECTOR3 v1(v0x,v0y,v0z);
				D3DXVECTOR3 v2(v2x,v2y,v2z);
				D3DXVECTOR3 v3(v1x,v1y,v1z);
				D3DXVECTOR3 vNorm;
				CD3DUtilityFunction::ComputeTriangleNormal(vNorm,v1,v2,v3);
				m_arrTriNormals.Add(vNorm);
			}
#endif
			//add triangle into the scene object			
			m_pScene->addGeometryObject(tri1);
			numTris++;
		}

		timer.Stop();
		double millsecs = timer.GetDuration();
		UtilityString strDebug;
		strDebug.GenerateFormatSmallString("%g milliseconds to createGeometryFrom2DVertArr()\n",millsecs);
		OutputDebugInfoThread(strDebug);

		return true;
}

bool CRayTracingHM::createRaysIn2DPlane(const D3DXVECTOR2& vLeftTopPt, float f2dDataWidth, float f2dDataHeight, float fInitZVal, float fScaleRatio,
	const D3DXVECTOR3& vRayDir)
{
	CHiPerfTimer timer;
	timer.Start();

	//compute how many rays per column
	unsigned int iNumRaysPerColIn2D = (unsigned int)ceil((f2dDataHeight/f2dDataWidth)*m_iNumRaysPerRowIn2D);

	D3DXVECTOR2 vNewLeftTopPt = vLeftTopPt + D3DXVECTOR2(0.5f*f2dDataWidth*(1.0f-fScaleRatio),0.5f*f2dDataHeight*(1.0f-fScaleRatio));
	float fNew2dDataWidth = f2dDataWidth * fScaleRatio;
	float fNew2dDataHeight = f2dDataHeight * fScaleRatio;
	float fStepX = (fNew2dDataWidth/m_iNumRaysPerRowIn2D);
	float fStepY = (fNew2dDataHeight/iNumRaysPerColIn2D);	

	//start to generate rays
	m_arrRays.RemoveAll();
	for (unsigned int j = 0; j < iNumRaysPerColIn2D; j++){
		float yPos = (j+0.5f)*fStepY + vNewLeftTopPt.y;
		for (unsigned int i = 0; i < m_iNumRaysPerRowIn2D; i++){
			float xPos = (i+0.5f)*fStepX + vNewLeftTopPt.x;

			Ray r;
			r.originX() = xPos; r.originY() = yPos; r.originZ() = fInitZVal;
			r.dirX() = vRayDir.x; r.dirY() = vRayDir.y; r.dirZ() = vRayDir.z;
			//from negative max to positive max
			r.tStart = -FLT_MAX; r.tEnd = FLT_MAX;

			m_arrRays.Add(r);
		}
	}
	timer.Stop();
	double millsecs = timer.GetDuration();
	UtilityString strDebug;
	strDebug.GenerateFormatSmallString("%g milliseconds to createRaysIn2DPlane()\n",millsecs);
	OutputDebugInfoThread(strDebug);

	return true;

}

bool CRayTracingHM::buildSceneHierarchy()
{
	if(!m_pScene)
		return false;

	CHiPerfTimer timer;
	timer.Start();
	switch(m_iCurHierBldStyle)
	{
	case FlatList:
		m_pScene->buildFlatList();
		break;
	case NaiveMiddleSplit:
		m_pScene->buildNaiveMiddleSplit(15);
		break;
	case NaiveMedianSplit:
		m_pScene->buildNaiveMedianSplit(8);
		break;
	case GreedySplit:
		m_pScene->buildGreedySplit();
		break;
	}
	timer.Stop();
	double millsecs = timer.GetDuration();

	UtilityString strDebug;
	strDebug.GenerateFormatSmallString("%g milliseconds to build acceleration structure\n",millsecs);
	OutputDebugInfoThread(strDebug);
	
	m_pScene->printTreeStats();

	return true;
}

bool CRayTracingHM::Tracing2dHeightMapData(const CXYZFileIO& curHeightMapFile, float fScaleRatio, const D3DXVECTOR3& vRayDir, const char* strOutputFN)
{
	unsigned int iWidth = curHeightMapFile.m_iHMDataWidth;
	unsigned int iHeight = curHeightMapFile.m_iHMDataHeight;
	//create geometry
	bool bRt = createGeometryFrom2DVertArr(curHeightMapFile.m_arrVertices, iWidth, iHeight);
	if(!bRt)
		return bRt;

	//build hierarchy
	bRt = buildSceneHierarchy();
	if(!bRt)
		return bRt;
	Triangle::printSlotStats();

	//create rays
	float fMaxHeight = curHeightMapFile.m_fHeightRange;
	float fInitialZVal = INITZVALSCALE*fMaxHeight;
	D3DXVECTOR3 vFirstPt = curHeightMapFile.m_arrVertices[0].Position;
	D3DXVECTOR2 vLeftTopPt(vFirstPt.x, vFirstPt.y);
	float f2dDataWidth = curHeightMapFile.m_fSampleWidth;
	float f2dDataHeight = curHeightMapFile.m_fSampleHeight;

	bRt = createRaysIn2DPlane(vLeftTopPt, f2dDataWidth, f2dDataHeight, fInitialZVal, fScaleRatio, vRayDir);
	if(!bRt)
		return bRt;
	int iTotalNumberOfRays = m_arrRays.GetSize();

	//start tracing with separate ray casting steps
	//for ping-pong ops
	CGrowableArray<Ray> arrRaysTmp;

	CGrowableArray<Ray>* pArrSrc = &m_arrRays;
	CGrowableArray<Ray>* pArrDst = &arrRaysTmp;

	static int iBounce = 0;
	FILE* fp = NULL;
	fopen_s(&fp,strOutputFN,"w");
	fprintf(fp, "TotalNumberOfRays=%d\n", iTotalNumberOfRays);
	fprintf(fp, "RayDir:(%f,%f,%f)\n", vRayDir.x, vRayDir.y, vRayDir.z);

	CHiPerfTimer timer;
	timer.Start();

	while(true)
	{
		//No Rays for tracing anymore
		if(pArrSrc->GetSize() == 0)
			break;

		fprintf(fp, "iBounce=%d, ", iBounce);
		//remove dst array
		pArrDst->RemoveAll();

		//number of rays that hit the 'No-data' points
		int iNumTracedRay = pArrSrc->GetSize();
		int iNumHitNoDataRay = 0;
		int iNumPointDownRay = 0;
		int iNumReflectOutRay = 0;
		int iNumRecordedRay = 0;
		int iNumMissedRay = 0;

		for (int iRay = 0; iRay < pArrSrc->GetSize(); iRay++)
		{
			IntersectionRecord hitRec;
			hitRec.primitiveId = -1;

			Ray ray = (*pArrSrc)[iRay];
			m_pScene->intersect(ray,&hitRec);

    	    if(hitRec.primitiveId == -1) {
				//no hit
				if(iBounce > 0)
					iNumReflectOutRay++;
				else
					iNumMissedRay++;				
			}
			else if(hitRec.primitiveId == 0)
				//'No-data' pt --> ignore for counting reflection
				iNumHitNoDataRay++;
			else{
				//find the normal of the hit triangle (id-1)!!!
				D3DXVECTOR3 vNorm = m_arrTriNormals[hitRec.primitiveId-1];
				//compute reflection vector
				D3DXVECTOR3 vIncident = D3DXVECTOR3(ray.dirX(),ray.dirY(),ray.dirZ());
				float vDotN = D3DXVec3Dot(&vNorm,&vIncident);
				D3DXVECTOR3 vTmp = vIncident - 2*vDotN*vNorm;
				D3DXVECTOR3 vReflect;
				D3DXVec3Normalize(&vReflect,&vTmp);

				//point-up reflection!
				if(vReflect.z > 0.0f) {
					//generate new ray
					D3DXVECTOR3 newOrgPt(ray.originX(), ray.originY(), ray.originZ());
					D3DXVECTOR3 vDir(ray.dirX(), ray.dirY(), ray.dirZ());
					newOrgPt += vDir * hitRec.t;

					Ray newray;
					newray.originX() = newOrgPt.x; newray.originY() = newOrgPt.y; newray.originZ() = newOrgPt.z;
					newray.dirX() = vReflect.x; newray.dirY() = vReflect.y; newray.dirZ() = vReflect.z;
					newray.tStart = ExtraEpsilon; newray.tEnd = TENDVAL - ExtraEpsilon;

					pArrDst->Add(newray);
				} else 
					iNumPointDownRay++;
			}
		}

		//output statistical data
		iNumRecordedRay = pArrDst->GetSize();
		fprintf(fp,"iNumTracedRay=%d, iNumReflectOutRay=%d, iNumPointDownRay=%d, iNumHitNoDataRay=%d, iNumRecordedRay=%d, iNumMissedRay=%d\n", 
			iNumTracedRay, iNumReflectOutRay, iNumPointDownRay, iNumHitNoDataRay, iNumRecordedRay, iNumMissedRay);

		//swap the two Src/Dst arrays
		CGrowableArray<Ray>* pTmp = pArrSrc;
		pArrSrc = pArrDst;
		pArrDst = pTmp;
		//increase bounce number
		iBounce++;
	}
	timer.Stop();
	double millsecs = timer.GetDuration();

	UtilityString strDebug;
	strDebug.GenerateFormatSmallString("%g milliseconds to finish ray tracing\n",millsecs);
	OutputDebugInfoThread(strDebug);

	//close file
	fclose(fp);
	iBounce = 0;
	
	return true;
}


bool CRayTracingHM::Tracing2dHeightMapData_RayOrder(const CXYZFileIO& curHeightMapFile, const float fScaleRatio, const D3DXVECTOR3& vRayDir, const char* strOutputFN)
{
	unsigned int iWidth = curHeightMapFile.m_iHMDataWidth;
	unsigned int iHeight = curHeightMapFile.m_iHMDataHeight;
	//create geometry
	bool bRt = createGeometryFrom2DVertArr(curHeightMapFile.m_arrVertices, iWidth, iHeight);
	if(!bRt)
		return bRt;

	//build hierarchy
	bRt = buildSceneHierarchy();
	if(!bRt)
		return bRt;
	Triangle::printSlotStats();

	//create rays
	float fMaxHeight = curHeightMapFile.m_fHeightRange;
	float fInitialZVal = INITZVALSCALE*fMaxHeight;
	D3DXVECTOR3 vFirstPt = curHeightMapFile.m_arrVertices[0].Position;
	D3DXVECTOR2 vLeftTopPt(vFirstPt.x, vFirstPt.y);
	float f2dDataWidth = curHeightMapFile.m_fSampleWidth;
	float f2dDataHeight = curHeightMapFile.m_fSampleHeight;

	bRt = createRaysIn2DPlane(vLeftTopPt, f2dDataWidth, f2dDataHeight, fInitialZVal, fScaleRatio, vRayDir);
	if(!bRt)
		return bRt;
	int iTotalNumberOfRays = m_arrRays.GetSize();

	//start tracing with separate ray casting steps
	FILE* fp = NULL;
	fopen_s(&fp,strOutputFN,"w");
	fprintf(fp, "TotalNumberOfRays=%d\n", iTotalNumberOfRays);
	fprintf(fp, "RayDir:(%f,%f,%f)\n", vRayDir.x, vRayDir.y, vRayDir.z);

	//for recording ray status
	CGrowableArray<RAYBOUNCETYPE> arrRayResults;
	CHiPerfTimer timer;
	timer.Start();

	for (int iRay = 0; iRay < m_arrRays.GetSize(); iRay++)
	{
		Ray ray = m_arrRays[iRay];
		int iBounce = 0;
		while (true)
		{
			IntersectionRecord hitRec;
			hitRec.primitiveId = -1;

			m_pScene->intersect(ray,&hitRec);

			if(hitRec.primitiveId == -1) {
				//no hit
				if(iBounce == 0)
					arrRayResults.Add(RAYMISSEDGEOMETRY);
				else if(iBounce == 1){
					if(ray.dirZ() > 0)
						arrRayResults.Add(RAYFIRSTBOUNCE);
					else
						arrRayResults.Add(RAYMISSEDGEOMETRY);
				}
				else if(iBounce > 1)
					arrRayResults.Add(RAYMULTIPLEBOUNCES);
				//quit loop
				break;
			} else if(hitRec.primitiveId == 0) {
				//'No-data' pt --> ignore for counting reflection
				if(iBounce == 0)
					arrRayResults.Add(RAYHITNODATAPT);
				else if(iBounce == 1)
					arrRayResults.Add(RAYHITNODATAPT_MBOUNCES);
				else 
					arrRayResults.Add(RAYMULTIPLEBOUNCES);
				//quit loop
				break;
			} else {
				//find the normal of the hit triangle (id-1)!!!
				D3DXVECTOR3 vNorm = m_arrTriNormals[hitRec.primitiveId-1];
				//compute reflection vector
				D3DXVECTOR3 vIncident = D3DXVECTOR3(ray.dirX(),ray.dirY(),ray.dirZ());
				float vDotN = D3DXVec3Dot(&vNorm,&vIncident);
				D3DXVECTOR3 vTmp = vIncident - 2*vDotN*vNorm;
				D3DXVECTOR3 vReflect;
				D3DXVec3Normalize(&vReflect,&vTmp);

				//generate new ray
				D3DXVECTOR3 newOrgPt(ray.originX(), ray.originY(), ray.originZ());
				D3DXVECTOR3 vDir(ray.dirX(), ray.dirY(), ray.dirZ());
				newOrgPt += vDir * hitRec.t;

				Ray newray;
				newray.originX() = newOrgPt.x; newray.originY() = newOrgPt.y; newray.originZ() = newOrgPt.z;
				newray.dirX() = vReflect.x; newray.dirY() = vReflect.y; newray.dirZ() = vReflect.z;
				newray.tStart = ExtraEpsilon; newray.tEnd = TENDVAL - ExtraEpsilon;
				//assign newray to ray!
				ray = newray;
			}
			//increase bounce
			iBounce++;
		}

	}

	timer.Stop();
	double millsecs = timer.GetDuration();

	UtilityString strDebug;
	strDebug.GenerateFormatSmallString("%g milliseconds to finished ray tracing\n",millsecs);
	OutputDebugInfoThread(strDebug);

	//count the number of rays for different statuses!
	int iNumMissedRay = 0;
	int iNumHitNoDataRay = 0;
	int iNumHitNoDataRayMbounces = 0;
	int iNumFirstBounceRef = 0;
	int iNumMultiBouncesRef = 0;
	for (int i = 0; i < arrRayResults.GetSize(); i++){
		switch(arrRayResults[i]){
		case RAYMISSEDGEOMETRY:
			iNumMissedRay++;
			break;
		case RAYHITNODATAPT:
			iNumHitNoDataRay++;
			break;
		case RAYHITNODATAPT_MBOUNCES:
			iNumHitNoDataRayMbounces++;
			break;
		case RAYFIRSTBOUNCE:
			iNumFirstBounceRef++;
			break;
		case RAYMULTIPLEBOUNCES:
			iNumMultiBouncesRef++;
			break;
		}
	}

	fprintf(fp,"iNumTracedRay=%d, iNumFirstBounceRef=%d, iNumMultiBouncesRef=%d, iNumHitNoDataRay=%d, iNumHitNoDataRayMbounces=%d, iNumMissedRay=%d\n", 
		iTotalNumberOfRays, iNumFirstBounceRef, iNumMultiBouncesRef, iNumHitNoDataRay, iNumHitNoDataRayMbounces, iNumMissedRay);
	fprintf(fp,"FirstBounceRatio=%f\n", (float)(iNumFirstBounceRef)/(iTotalNumberOfRays-iNumHitNoDataRay));
	//close file
	fclose(fp);

	strDebug.GenerateFormatSmallString("iNumTracedRay=%d, iNumFirstBounceRef=%d, iNumMultiBouncesRef=%d, iNumHitNoDataRay=%d, iNumHitNoDataRayMbounces=%d, iNumMissedRay=%d\n", 
		iTotalNumberOfRays, iNumFirstBounceRef, iNumMultiBouncesRef, iNumHitNoDataRay, iNumHitNoDataRayMbounces, iNumMissedRay);
	OutputDebugInfoThread(strDebug);
	strDebug.GenerateFormatSmallString("FirstBounceRatio=%f\n", (float)(iNumFirstBounceRef)/(iTotalNumberOfRays-iNumHitNoDataRay));
	OutputDebugInfoThread(strDebug);

	return true;
}


bool CRayTracingHM::Tracing2dHMDataVirutalGonio_OneOfMultiXYZFiles(const CXYZFileIO& curHeightMapFile, const float fScaleRatio, const D3DXVECTOR3& vRayDir, 
	double* arrSumBuf, double* arrSingleBuf, double* arrMultiBuf, unsigned int iArrBufWidth, unsigned int iArrBufHeight,
	const double iorIncidentReal, const double iorTransmittedReal, const double iorTransmittedImaginary)
{
	if(!arrSumBuf || !arrSingleBuf || !arrMultiBuf)
		return false;

	unsigned int iWidth = curHeightMapFile.m_iHMDataWidth;
	unsigned int iHeight = curHeightMapFile.m_iHMDataHeight;
	//create geometry
	bool bRt = createGeometryFrom2DVertArr(curHeightMapFile.m_arrVertices, iWidth, iHeight);
	if(!bRt)
		return bRt;

	//build hierarchy
	bRt = buildSceneHierarchy();
	if(!bRt)
		return bRt;
	//Triangle::printSlotStats();

	//create rays
	float fMaxHeight = curHeightMapFile.m_fHeightRange;
	float fInitialZVal = INITZVALSCALE*fMaxHeight;
	D3DXVECTOR3 vFirstPt = curHeightMapFile.m_arrVertices[0].Position;
	D3DXVECTOR2 vLeftTopPt(vFirstPt.x, vFirstPt.y);
	float f2dDataWidth = curHeightMapFile.m_fSampleWidth;
	float f2dDataHeight = curHeightMapFile.m_fSampleHeight;

	bRt = createRaysIn2DPlane(vLeftTopPt, f2dDataWidth, f2dDataHeight, fInitialZVal, fScaleRatio, vRayDir);
	if(!bRt)
		return bRt;
	int iTotalNumberOfRays = m_arrRays.GetSize();
	//the energy for each ray
	double fEnergyEachRay = m_fTotalIncidentEnergy/iTotalNumberOfRays;

	//start tracing with separate ray casting steps
	//for recording ray status
	CGrowableArray<RAYBOUNCETYPE> arrRayResults;
	//array for saving the exitant energy
	CHiPerfTimer timer;
	timer.Start();	

	for (int iRay = 0; iRay < m_arrRays.GetSize(); iRay++)
	{
		Ray ray = m_arrRays[iRay];
		int iBounce = 0;
		double fEnergyCurRay = fEnergyEachRay;
		while (true)
		{
			IntersectionRecord hitRec;
			hitRec.primitiveId = -1;

			m_pScene->intersect(ray,&hitRec);

			if(hitRec.primitiveId == -1) {
				//no hit
				if(iBounce == 0)
					arrRayResults.Add(RAYMISSEDGEOMETRY);
				else if(iBounce == 1){
					if(ray.dirZ() > 0) {
						arrRayResults.Add(RAYFIRSTBOUNCE);
						D3DXVECTOR3 vVec = D3DXVECTOR3(ray.dirX(),ray.dirY(),ray.dirZ());
						//save energy
						SaveTo2dGonioBuf(vVec,fEnergyCurRay,arrSumBuf,iArrBufWidth,iArrBufHeight);
						SaveTo2dGonioBuf(vVec,fEnergyCurRay,arrSingleBuf,iArrBufWidth,iArrBufHeight);
					}else
						arrRayResults.Add(RAYMISSEDGEOMETRY);

				} else if(iBounce > 1){
					if(ray.dirZ() > 0){
						arrRayResults.Add(RAYMULTIPLEBOUNCES);
						D3DXVECTOR3 vVec = D3DXVECTOR3(ray.dirX(),ray.dirY(),ray.dirZ());
						//save energy
						SaveTo2dGonioBuf(vVec,fEnergyCurRay,arrSumBuf,iArrBufWidth,iArrBufHeight);
						SaveTo2dGonioBuf(vVec,fEnergyCurRay,arrMultiBuf,iArrBufWidth,iArrBufHeight);
					} else 
						arrRayResults.Add(RAYMISSEDGEOMETRY);
				}
				//quit loop
				break;
			} else if(hitRec.primitiveId == 0){ 
				//'No-data' pt --> ignore for counting reflection
				if(iBounce == 0)
					arrRayResults.Add(RAYHITNODATAPT);
				else if(iBounce == 1)
					arrRayResults.Add(RAYHITNODATAPT_MBOUNCES);
				else 
					arrRayResults.Add(RAYMULTIPLEBOUNCES);
				//quit loop
				break;
			} else {
				//find the normal of the hit triangle (id-1)!!!
				D3DXVECTOR3 vNorm = m_arrTriNormals[hitRec.primitiveId-1];
				//compute reflection vector
				D3DXVECTOR3 vIncident = D3DXVECTOR3(ray.dirX(),ray.dirY(),ray.dirZ());
				float vDotN = D3DXVec3Dot(&vNorm,&vIncident);
				D3DXVECTOR3 vTmp = vIncident - 2*vDotN*vNorm;
				D3DXVECTOR3 vReflect;
				D3DXVec3Normalize(&vReflect,&vTmp);

				//generate new ray
				D3DXVECTOR3 newOrgPt(ray.originX(), ray.originY(), ray.originZ());
				D3DXVECTOR3 vDir(ray.dirX(), ray.dirY(), ray.dirZ());
				newOrgPt += vDir * hitRec.t;

				Ray newray;
				newray.originX() = newOrgPt.x; newray.originY() = newOrgPt.y; newray.originZ() = newOrgPt.z;
				newray.dirX() = vReflect.x; newray.dirY() = vReflect.y; newray.dirZ() = vReflect.z;
				newray.tStart = ExtraEpsilon; newray.tEnd = TENDVAL - ExtraEpsilon;
				//assign newray to ray!
				ray = newray;

				//update the energy for the reflected ray!
				vDotN *= -1.0;//inverse cosine
				double fFresenel = CD3DUtilityFunction::computeFresnelComplex(iorIncidentReal,iorTransmittedReal,iorTransmittedImaginary,vDotN);
				fEnergyCurRay *= fFresenel;
				//fEnergyCurRay *= 1.0;
			}
			//increase bounce
			iBounce++;
		}

	}

	timer.Stop();
	double millsecs = timer.GetDuration();

	UtilityString strDebug;
	strDebug.GenerateFormatSmallString("%g milliseconds to finished ray tracing\n",millsecs);
	OutputDebugInfoThread(strDebug);

	//count the number of rays for different statuses!
	int iNumMissedRay = 0;
	int iNumHitNoDataRay = 0;
	int iNumHitNoDataRayMbounces = 0;
	int iNumFirstBounceRef = 0;
	int iNumMultiBouncesRef = 0;
	for (int i = 0; i < arrRayResults.GetSize(); i++){
		switch(arrRayResults[i]){
		case RAYMISSEDGEOMETRY:
			iNumMissedRay++;
			break;
		case RAYHITNODATAPT:
			iNumHitNoDataRay++;
			break;
		case RAYHITNODATAPT_MBOUNCES:
			iNumHitNoDataRayMbounces++;
			break;
		case RAYFIRSTBOUNCE:
			iNumFirstBounceRef++;
			break;
		case RAYMULTIPLEBOUNCES:
			iNumMultiBouncesRef++;
			break;
		}
	}

	strDebug.GenerateFormatSmallString("iNumTracedRay=%d, iNumFirstBounceRef=%d, iNumMultiBouncesRef=%d, iNumHitNoDataRay=%d, iNumHitNoDataRayMbounces=%d, iNumMissedRay=%d\n", 
		iTotalNumberOfRays, iNumFirstBounceRef, iNumMultiBouncesRef, iNumHitNoDataRay, iNumHitNoDataRayMbounces, iNumMissedRay);
	OutputDebugInfoThread(strDebug);
	strDebug.GenerateFormatSmallString("FirstBounceRatio=%f\n", (float)(iNumFirstBounceRef)/(iTotalNumberOfRays-iNumHitNoDataRay));
	OutputDebugInfoThread(strDebug);
	
	return true;
}



bool CRayTracingHM::Tracing2dHMDataVirutalGonio(const CXYZFileIO& curHeightMapFile, const float fScaleRatio, const D3DXVECTOR3& vRayDir, 
	const char* strOutputFN, const char* strSumFN, const char* strFirstFN, const char* strMultiFN, 
	const double iorIncidentReal, const double iorTransmittedReal, const double iorTransmittedImaginary)
{
	unsigned int iWidth = curHeightMapFile.m_iHMDataWidth;
	unsigned int iHeight = curHeightMapFile.m_iHMDataHeight;
	//create geometry
	bool bRt = createGeometryFrom2DVertArr(curHeightMapFile.m_arrVertices, iWidth, iHeight);
	if(!bRt)
		return bRt;

	//build hierarchy
	bRt = buildSceneHierarchy();
	if(!bRt)
		return bRt;
	//Triangle::printSlotStats();

	//create rays
	float fMaxHeight = curHeightMapFile.m_fHeightRange;
	float fInitialZVal = INITZVALSCALE*fMaxHeight;
	D3DXVECTOR3 vFirstPt = curHeightMapFile.m_arrVertices[0].Position;
	D3DXVECTOR2 vLeftTopPt(vFirstPt.x, vFirstPt.y);
	float f2dDataWidth = curHeightMapFile.m_fSampleWidth;
	float f2dDataHeight = curHeightMapFile.m_fSampleHeight;

	bRt = createRaysIn2DPlane(vLeftTopPt, f2dDataWidth, f2dDataHeight, fInitialZVal, fScaleRatio, vRayDir);
	if(!bRt)
		return bRt;
	int iTotalNumberOfRays = m_arrRays.GetSize();
	//the energy for each ray
	double fEnergyEachRay = m_fTotalIncidentEnergy/iTotalNumberOfRays;
	
	//start tracing with separate ray casting steps
	FILE* fp = NULL;
	fopen_s(&fp,strOutputFN,"w");
	fprintf(fp, "TotalNumberOfRays=%d\n", iTotalNumberOfRays);
	fprintf(fp, "RayDir:(%f,%f,%f)\n", vRayDir.x, vRayDir.y, vRayDir.z);

	//for recording ray status
	CGrowableArray<RAYBOUNCETYPE> arrRayResults;
	//array for saving the exitant energy
	double* arrExitEnergy = new double[m_iVGonioOutputImgWidth*m_iVGonioOutputImgHeight];
	memset(arrExitEnergy,0,sizeof(double)*m_iVGonioOutputImgWidth*m_iVGonioOutputImgHeight);
	double* arrExitEnergy_Single = new double[m_iVGonioOutputImgWidth*m_iVGonioOutputImgHeight];
	memset(arrExitEnergy_Single,0,sizeof(double)*m_iVGonioOutputImgWidth*m_iVGonioOutputImgHeight);
	double* arrExitEnergy_Multiple = new double[m_iVGonioOutputImgWidth*m_iVGonioOutputImgHeight];
	memset(arrExitEnergy_Multiple,0,sizeof(double)*m_iVGonioOutputImgWidth*m_iVGonioOutputImgHeight);

	CHiPerfTimer timer;
	timer.Start();	

	for (int iRay = 0; iRay < m_arrRays.GetSize(); iRay++)
	{
		Ray ray = m_arrRays[iRay];
		int iBounce = 0;
		double fEnergyCurRay = fEnergyEachRay;
		while (true)
		{
			IntersectionRecord hitRec;
			hitRec.primitiveId = -1;

			m_pScene->intersect(ray,&hitRec);

			if(hitRec.primitiveId == -1) {
				//no hit
				if(iBounce == 0)
					arrRayResults.Add(RAYMISSEDGEOMETRY);
				else if(iBounce == 1){
					if(ray.dirZ() > 0) {
						arrRayResults.Add(RAYFIRSTBOUNCE);
						D3DXVECTOR3 vVec = D3DXVECTOR3(ray.dirX(),ray.dirY(),ray.dirZ());
						//save energy
						SaveTo2dGonioBuf(vVec,fEnergyCurRay,arrExitEnergy,m_iVGonioOutputImgWidth,m_iVGonioOutputImgHeight);
						SaveTo2dGonioBuf(vVec,fEnergyCurRay,arrExitEnergy_Single,m_iVGonioOutputImgWidth,m_iVGonioOutputImgHeight);
					}else
						arrRayResults.Add(RAYMISSEDGEOMETRY);
					
				} else if(iBounce > 1){
					if(ray.dirZ() > 0){
						arrRayResults.Add(RAYMULTIPLEBOUNCES);
						D3DXVECTOR3 vVec = D3DXVECTOR3(ray.dirX(),ray.dirY(),ray.dirZ());
						//save energy
						SaveTo2dGonioBuf(vVec,fEnergyCurRay,arrExitEnergy,m_iVGonioOutputImgWidth,m_iVGonioOutputImgHeight);
						SaveTo2dGonioBuf(vVec,fEnergyCurRay,arrExitEnergy_Multiple,m_iVGonioOutputImgWidth,m_iVGonioOutputImgHeight);
					} else 
						arrRayResults.Add(RAYMISSEDGEOMETRY);
				}
				//quit loop
				break;
			} else if(hitRec.primitiveId == 0){ 
				//'No-data' pt --> ignore for counting reflection
				if(iBounce == 0)
					arrRayResults.Add(RAYHITNODATAPT);
				else if(iBounce == 1)
					arrRayResults.Add(RAYHITNODATAPT_MBOUNCES);
				else 
					arrRayResults.Add(RAYMULTIPLEBOUNCES);
				//quit loop
				break;
			} else {
				//find the normal of the hit triangle (id-1)!!!
				D3DXVECTOR3 vNorm = m_arrTriNormals[hitRec.primitiveId-1];
				//compute reflection vector
				D3DXVECTOR3 vIncident = D3DXVECTOR3(ray.dirX(),ray.dirY(),ray.dirZ());
				float vDotN = D3DXVec3Dot(&vNorm,&vIncident);
				D3DXVECTOR3 vTmp = vIncident - 2*vDotN*vNorm;
				D3DXVECTOR3 vReflect;
				D3DXVec3Normalize(&vReflect,&vTmp);

				//generate new ray
				D3DXVECTOR3 newOrgPt(ray.originX(), ray.originY(), ray.originZ());
				D3DXVECTOR3 vDir(ray.dirX(), ray.dirY(), ray.dirZ());
				newOrgPt += vDir * hitRec.t;

				Ray newray;
				newray.originX() = newOrgPt.x; newray.originY() = newOrgPt.y; newray.originZ() = newOrgPt.z;
				newray.dirX() = vReflect.x; newray.dirY() = vReflect.y; newray.dirZ() = vReflect.z;
				newray.tStart = ExtraEpsilon; newray.tEnd = TENDVAL - ExtraEpsilon;
				//assign newray to ray!
				ray = newray;

				//update the energy for the reflected ray!
				vDotN *= -1.0;//inverse cosine
				double fFresenel = CD3DUtilityFunction::computeFresnelComplex(iorIncidentReal,iorTransmittedReal,iorTransmittedImaginary,vDotN);
				//fEnergyCurRay *= fFresenel;
				fEnergyCurRay *= 1.0;
			}
			//increase bounce
			iBounce++;
		}

	}

	timer.Stop();
	double millsecs = timer.GetDuration();

	UtilityString strDebug;
	strDebug.GenerateFormatSmallString("%g milliseconds to finished ray tracing\n",millsecs);
	OutputDebugInfoThread(strDebug);

	//count the number of rays for different statuses!
	int iNumMissedRay = 0;
	int iNumHitNoDataRay = 0;
	int iNumHitNoDataRayMbounces = 0;
	int iNumFirstBounceRef = 0;
	int iNumMultiBouncesRef = 0;
	for (int i = 0; i < arrRayResults.GetSize(); i++){
		switch(arrRayResults[i]){
		case RAYMISSEDGEOMETRY:
			iNumMissedRay++;
			break;
		case RAYHITNODATAPT:
			iNumHitNoDataRay++;
			break;
		case RAYHITNODATAPT_MBOUNCES:
			iNumHitNoDataRayMbounces++;
			break;
		case RAYFIRSTBOUNCE:
			iNumFirstBounceRef++;
			break;
		case RAYMULTIPLEBOUNCES:
			iNumMultiBouncesRef++;
			break;
		}
	}

	fprintf(fp,"iNumTracedRay=%d, iNumFirstBounceRef=%d, iNumMultiBouncesRef=%d, iNumHitNoDataRay=%d, iNumHitNoDataRayMbounces=%d, iNumMissedRay=%d\n", 
		iTotalNumberOfRays, iNumFirstBounceRef, iNumMultiBouncesRef, iNumHitNoDataRay, iNumHitNoDataRayMbounces, iNumMissedRay);
	fprintf(fp,"FirstBounceRatio=%f\n", (float)(iNumFirstBounceRef)/(iTotalNumberOfRays-iNumHitNoDataRay));
	//close file
	fclose(fp);

	strDebug.GenerateFormatSmallString("iNumTracedRay=%d, iNumFirstBounceRef=%d, iNumMultiBouncesRef=%d, iNumHitNoDataRay=%d, iNumHitNoDataRayMbounces=%d, iNumMissedRay=%d\n", 
		iTotalNumberOfRays, iNumFirstBounceRef, iNumMultiBouncesRef, iNumHitNoDataRay, iNumHitNoDataRayMbounces, iNumMissedRay);
	OutputDebugInfoThread(strDebug);
	strDebug.GenerateFormatSmallString("FirstBounceRatio=%f\n", (float)(iNumFirstBounceRef)/(iTotalNumberOfRays-iNumHitNoDataRay));
	OutputDebugInfoThread(strDebug);

	//save virtual gonio output results!
	//find the maximum value
	double fMaxVal = -FLT_MAX;
	for (unsigned int j = 0; j < m_iVGonioOutputImgHeight; j++)
		for (unsigned int i = 0; i < m_iVGonioOutputImgWidth; i++)
		{
			double val = arrExitEnergy[j*m_iVGonioOutputImgWidth+i];
			if(val > fMaxVal)
				fMaxVal = val;
		}
#if 0	
	Save2dGonioBufTo8BitImg(arrExitEnergy,m_iVGonioOutputImgWidth,m_iVGonioOutputImgHeight,strSumFN,0,fMaxVal,false);
	strDebug.GenerateFormatSmallString("%s saved!", strSumFN);
	OutputDebugInfoThread(strDebug);
	Save2dGonioBufTo8BitImg(arrExitEnergy_Single,m_iVGonioOutputImgWidth,m_iVGonioOutputImgHeight,strFirstFN,0,fMaxVal,false);
	strDebug.GenerateFormatSmallString("%s saved!", strFirstFN);
	OutputDebugInfoThread(strDebug);
	Save2dGonioBufTo8BitImg(arrExitEnergy_Multiple,m_iVGonioOutputImgWidth,m_iVGonioOutputImgHeight,strMultiFN,0,fMaxVal,false);
	strDebug.GenerateFormatSmallString("%s saved!", strMultiFN);
	OutputDebugInfoThread(strDebug);
#else
	Save2dGonioBufTo16BitImg(arrExitEnergy,m_iVGonioOutputImgWidth,m_iVGonioOutputImgHeight,strSumFN,0,fMaxVal,false);
	strDebug.GenerateFormatSmallString("%s saved!", strSumFN);
	OutputDebugInfoThread(strDebug);
	Save2dGonioBufTo16BitImg(arrExitEnergy_Single,m_iVGonioOutputImgWidth,m_iVGonioOutputImgHeight,strFirstFN,0,fMaxVal,false);
	strDebug.GenerateFormatSmallString("%s saved!", strFirstFN);
	OutputDebugInfoThread(strDebug);
	Save2dGonioBufTo16BitImg(arrExitEnergy_Multiple,m_iVGonioOutputImgWidth,m_iVGonioOutputImgHeight,strMultiFN,0,fMaxVal,false);
	strDebug.GenerateFormatSmallString("%s saved!", strMultiFN);
	OutputDebugInfoThread(strDebug);
#endif
	//release array memory
	SAFE_DELETE_ARRAY(arrExitEnergy);
	SAFE_DELETE_ARRAY(arrExitEnergy_Single);
	SAFE_DELETE_ARRAY(arrExitEnergy_Multiple);

	return true;
}

void CRayTracingHM::SaveTo2dGonioBuf(const D3DXVECTOR3& vVec, const double fEnergyCurVec, double* arrStoreBuf, unsigned int iBufWidth, unsigned int iBufHeight)
{
	if(!arrStoreBuf)
		return;

	D3DXVECTOR3 vec;
	D3DXVec3Normalize(&vec, &vVec);
	vec.x = vec.x/2.0f + 0.5f; //(-1,1) --> (0,1); 
	vec.y = vec.y/2.0f + 0.5f;
	int iPosX = (int)(floorf(vec.x * iBufWidth));
	int iPosY = (int)(floorf(vec.y * iBufHeight));
	if(iPosX == iBufWidth) iPosX--;
	if(iPosY == iBufHeight) iPosY--;
	unsigned int idx = iPosY*iBufWidth + iPosX;
	arrStoreBuf[idx] += fEnergyCurVec;
}

void CRayTracingHM::Save2dGonioBufTo8BitImg(const double* arrStoreBuf, unsigned int iBufWidth, unsigned int iBufHeight, const char* strImgFN, 
	const double fMinVal, const double fMaxVal, bool bColor)
{
	if(!arrStoreBuf)
		return;

	//prepare CImg object
		if(bColor){
		CImg<unsigned char> oImg(iBufWidth,iBufHeight,1,3,0);
		for (unsigned int j = 0; j < iBufHeight; j++)
			for (unsigned int i = 0; i < iBufWidth; i++)
			{
				double val = arrStoreBuf[j*iBufWidth+i];
				double fNormalizedVal = val/(fMaxVal-fMinVal);//[0,1]

				D3DXVECTOR3 color = CD3DUtilityFunction::GetColour(static_cast<float>(fNormalizedVal),0.0f,1.0f);
				oImg.atXY(i,j,0,0) = static_cast<unsigned char>(255*color.x);			
				oImg.atXY(i,j,0,1) = static_cast<unsigned char>(255*color.y);			
				oImg.atXY(i,j,0,2) = static_cast<unsigned char>(255*color.z);	
			}
		//save image	
		oImg.save_bmp(strImgFN);
		oImg.clear();
		}else {
			CImg<unsigned char> oImg(iBufWidth,iBufHeight,1,1,0);
			for (unsigned int j = 0; j < iBufHeight; j++)
				for (unsigned int i = 0; i < iBufWidth; i++)
				{
					double val = arrStoreBuf[j*iBufWidth+i];
					double fNormalizedVal = val/(fMaxVal-fMinVal);//[0,1]
					oImg.atXY(i,j,0,0) = static_cast<unsigned char>(255*fNormalizedVal);								
				}
			//save image	
			oImg.save_bmp(strImgFN);
			oImg.clear();
		}
			
}

void CRayTracingHM::Save2dGonioBufTo16BitImg(const double* arrStoreBuf, unsigned int iBufWidth, unsigned int iBufHeight, 
	const char* strImgFN, const double fMinVal, const double fMaxVal, bool bColor/* =true */)
{
	CExrFileIO* pExrFileIO = new CExrFileIO();

	pExrFileIO->FillPixelsArray(arrStoreBuf, iBufWidth, iBufHeight, (fMaxVal - fMinVal), false );
	pExrFileIO->SaveExrFile(const_cast<char*>(strImgFN), Imf::WRITE_Y);

	SAFE_DELETE(pExrFileIO);
}